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Abstract. We provide a comprehensive study on the planar (2D) orienta¬ 
tional distributions of nematic polymers under an imposed shear flow of ar¬ 
bitrary strength. We extend previous analysis for persistence of equilibria in 
steady shear and for transitions to unsteady limit cycles, from closure models 
[21] to the Doi-Hess 2D kinetic equation. A variation on the Boltzmann distri¬ 
bution analysis of Constantin et al. [3, 4, 5] and others [8, 22, 23] for potential 
flow is developed to solve for all persistent steady equilibria, and characterize 
parameter boundaries where steady states cease to exist, which predicts the 
transition to tumbling limit cycles. 


1. Introduction. Macromolecular materials (for example, nematic liquid crys¬ 
tal polymers (LCP)) exhibit large flexibility and can be easily processed into fibers 
with high strength, or thin films [1, 6, 26]. The bulk properties of these materials 
are closely related to the processing flows. So it is very important to understand 
the dynamic behavior of the materials in the presence of flows. 

The kinetic Doi-Hess theory [7, 17] has been a popular model for nematic LCPs. 
In kinetic theory each nematogenic molecule is coarse grained as a rigid rod and the 
ensemble is described by an orientational probability density function (pdf) which 
evolves according to the Smoluchowski (or Fokker-Planck) equation. Theoretical 
investigations have been carried out for pure nematic equilibria [3, 4, 5, 8, 22, 23, 
29], extensional flow-induced equilibria [27], stable equilibria of dipolar ensembles 
[19, 30], effect of high [24, 25, 28] and weak shear [31], and effect of coplanar magnetic 
field [16]. Numerical simulations of the Smoluchowski equations can be found in 
[2, 9, 10, 11, 12, 13, 14, 15, 18, 20]. The first 3-D simulations carried out by Larson 
[20] showed shear-induced transitions and confirmed Marrucci and Maffettone’s 2-D 
calculations [24, 25]. 
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In the 2-D case, it is well-known that in the absence of flow the isotropic-nematic 
phase transition occurs at U = 2, where U is the normalized polymer concentration. 
In the presence of an imposed weak shear there is a threshold (Uo ~ 2.41144646) for 
U: When U < Uq, steady state solution exists; otherwise there is no steady state 
and the orientational pdf is temporally periodic (“tumbling”) and can never reach 
a steady state. 

The goal of this work is to conduct a comprehensive study of the 2-D nematic 
liquid crystal polymers under an arbitrary shear. In [21] monolayer films of liquid 
crystal polymers have been modeled with a mesoscopic two-dimensional (2D) ana¬ 
logue of the Doi-Hess kinetic model. The weak-shear steady and unsteady selection 
criteria for 2D nematic polymers using various second-moment closures was de¬ 
rived. A simple proof was given based on the Poincare-Bendixon Theorem to show 
that limit cycles (“tumbling orbits”) exist beyond the parameter boundary for the 
steady-unsteady transition. Finally, it was revealed that the shear-perturbed 2D 
phase diagram is significantly robust to closure approximations than the 3D sys¬ 
tem. Our approach here differs from [21] in that we use the kinetic Smoluchowski 
equation instead of the Doi tensor model. There are two minor differences between 
our work and [21]. 1) In [21], a bistable region (where there are two stable solu¬ 
tions) was observed for the Doi closure tensor model. In our study using the full 
Smoluchowski equation there is no bistable region. Note that the bistable region 
in [21] is extremely small and it is difficult to pinpoint what small differences be¬ 
tween the full Smoluchowski model and closure-based tensor models might cause 
this small bistable region. 2) In [21], Bogdanov-Takens bifurcation was observed 
whereas there is no Bogdanov-Takens bifurcation in our study. Again note that in 
[21] the Bogdanov-Takens bifurcation occurs in a very narrow region in the phase 
diagram. Our results are in general consistent with those in [21] in several aspects: 
1) Both models identify critical values of polymer concentration, beyond which there 
is no steady state; 2) Both models yield the fold structure in the phase diagram. 
However, a detailed fold structure is resolved in this work. Specifically, the fold 
is not a kink, rather it is a smooth turn with large curvature. Furthermore, we 
find that fold structure in the phase diagram (the plot of order parameter versus 
U) exits for Pe < 0.746. For Pe > 0.746, the fold structure disappears. At the 
critical value of Pe « 0.746, the boundary separating the steady state region and 
the tumbing region in the Pe — U plane appears to have a singularity. In [21] it 
was mentioned that unstable state can be further classified as tumbling or wagging. 
Tensor models were originally motivated by the fact that in the Doi-Hess kinetic 
model with the Maier-Saupe potential, an equilibrium state is completely specified 
by the orientation tensor. As a result, tensor models can provide good approxima¬ 
tions at equilibrium or near equilibrium. In 2-D tensor models the system state 
is represented by the orientation direction and the order parameter. Thus, it is 
natural to try to classify time-periodic unsteady state as tumbling or wagging. If 
the director rotates continuously in one direction, we may say it is tumbling; if the 
director moves back and forth in a certain range, we may say it is wagging. Even in 
the tensor models, such a classification may be problematic since the director angle 
is undefined at the isotropic state. If the time-periodic evolution ever goes through 
the isotropic state, such a classification may be vague. It is worthwhile to note that 
“tumbling” refers to the tumbling of the orientation probability density. As for 
individual polymer rods, under a shear they are always tumbling in the direction 
of shear even if the orientation pdf is in a steady state. With a full Smoluchowski 
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equation (i.e. the Doi-Hess model),we find that under a weak shear the orientation 
probability density function (pdf) is a traveling wave with the location-dependent 
velocity, that is, the shape of the orientation pdf is maintained while the whole 
function may shift in the orientation dimension. So in the case of weak shear, 
the classification of tumbling and wagging makes perfect sense. But in the case 
of weak shear, it was found that there is only tumbling and there is no wagging 
[31]. As Pe gets larger (i.e. the shear getting stronger), the unsteady evolution of 
the orientation pdf becomes more complicated. As the peak of the pdf moves, the 
shape of the pdf changes dramatically for large Pe. We feel that the classification 
of tumbling and wagging is no longer adequate for providing a good picture of the 
time-evolution of the orientation pdf. We will pursue a comprehensive classification 
of the pdf in the future work. Here we shall refer the time-periodic unsteady state 
simply as tumbling. 

In this paper we present and discuss phase diagrams for a very wide range of 
the Peclet number. For small Peclet number, the result of the current analyt¬ 
ical/numerical study confirms the conclusions of the previous asymptotic study 
[31]: When the concentration parameter U is greater than a critical value Uq 
(« 2.41144646), there is no stable steady state; when U is smaller than Uq 7 there is 
one stable branch of steady state; the stable branch and an unstable branch form a 
fold and are connected at Uq. The current study resolves the details of the transition 
from stability to instability along the phase diagram: the transition occurs around 
the folding point of the phase curve. As the Peclet number increases, the stable 
branch and the unstable branch of the fold are peeled off from each other and are 
separated. For Pe > 0.746, the fold disappears completely, that is, the phase curve 
no longer has a folding point. Furthermore, for Pe > 0.746, there is only one steady 
state (stable or unstable) for each value of U, and the transition from stability to 
instability occurs at a Hopf bifurcation point. 

2. Two-Dimensional formulation. We start with reviewing the mathemat¬ 
ical formulation of the Doi-Hess kinetic theory for two-dimensional homogeneous 
flows of rigid rodlike molecules immersed in a viscous solvent subject to an imposed 
flow field [7, 17]. The orientation of a polymer rod is described by an angle 9 and 
the orientational direction of each polymer rod is denoted by u = (cos 6, sin 9). 

We consider the steady state solution of the Smoluchowski equation in the form 


[7] 



{-f + ip\9))p+ -^p 


(1) 


where p(9,t) is the orientational probability density function of the ensemble, ip (9) 
is a periodic function with period 7r, and / is a constant torque. In the case of 
nematic polymers under shear, ip{9) contains both the Maier-Saupe interaction and 
the effect of the elongational component of the shear whereas / is caused by the 
rotational component of the shear. The steady state solution satisfies 



( 2 ) 


where J is the steady state probability flux. Multiplying (2) by the integrating 
factor exp [—f9 + ip {6)], we have 


[exp(-/0 + ip(9))p\ = - J exp [-/ 9 + ip{9)}. 


(3) 
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Integrating from 0 to 0 + n yields 


[exp(—/ 7 r) - 1] exp[-/0 + ip(9)]p{6) = - J 


(•O + 'K 


exp [—f0 + ij}{9)\d9. 


(4) 


Solving for p{9) gives 

p(o) = 1 _ ex ^_ exp[-V , (^)] J exp [-/e+ ip{e+ e))de. (5) 

The value of the probability flux J is determined from the condition 

p{9)d9 = 1. ( 6 ) 

Note that once we know / and ip{9), the probability density p{9) is conveniently 
determined from (5). An efficient method of implementing (5) based on the Fast 
Fourier Transform (FFT) is discussed in Appendix A. 

When an external shear flow with arbitrary strength is imposed, the Smolu- 
chowski equation ( 1 ) can be written as 



dp 

dt 


d_ 

dd 


( y Vsh(@) + Vms(&))p + 


(7) 


where Pe = 'y/Dr is the Peclet number, a nondimensional parameter measuring 
the relative strength of viscosity ( 7 ) and rotational diffusivity {Dr). In equation 
(7), Vsh{0) is the periodic part of the potential induced by the imposed shear, 

V S h{0) = ^ sin 20 , ( 8 ) 

and Vms{9) is the Maier-Saupe short-range interaction potential, 


Vms{6) = —U(cos2{9 — a)) cos 2 {9 — a), (9) 


where the brackets (•) denote the ensemble average with respect to the probability 
density function p{9) and the angle a is selected such that 

(sin 2 ( 0 - a)) = 0 . ( 10 ) 


Let r = [/(cos 2(0 — a)). With the introduction of r, the Maier-Saupe potential can 
be expressed as 

Vms {6) = -r cos 2(0 — a) ( 11 ) 

and the order parameter can be written as 

(cos 2(0 — a)) = jj. ( 12 ) 

Note that equation ( 1 ) reduces to (7) if we use the substitutions / = ^ and 


m 


Pe 

— sin 20 — r cos 2 (9 — a) 

4 

Pc Pc 

— cos 2a sin 2(0 — a) — (r-— sin 2a) cos 2(0 — a) 

Pe 

— cos 2a sin 2(0 — a) — geos 2(0 — a) 


(13) 


where 


q = r — 


„ 

— sm 2a. 
4 


(14) 


It is clear that once the values of Pe, r and a are given, the total potential is 
completely specified and the steady state probability density (if a steady state exists) 
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is given by (5). The existence of steady state solution of the Smoluchowski equation 
(7) depends on the existence of solution of the following nonlinear system: 

Pl(o, q,Pe) = (sin 2(6 — a)) = 0, 

F 2 (a,q,Pe) = (cos 2 (6 - a)) = (15) 

where the pdf used in averaging is the one given in (5). In the nonlinear system 
(15), Pe and U are parameters whereas r and a are unknowns. Alternatively, q and 
a can be used as unknowns where q and r are related by (14). 

The goal of this paper is to study the existence and stability of steady solutions 
of the nonlinear system (15) for Pe > 0 and U > 0. It is worth noticing that with 
the introduction of r, the pdf given in (5) is independent of U. So we can avoid the 
parameter U in our analysis until the very end. These are the big advantage and 
mathematical convenience of dealing with r instead of U. 


3. Existence of steady state solutions. In this section we will study (both 
analytically and numerically) the existence of solutions of the nonlinear system (15) 
with 0 < U < oo for a given value of Pe. It is important to keep in mind that the 
pdf in (5) is independent of U with the use of r. 

Since sin 2 6 is a periodic function with period ir, we can restrict the range of a 
to [— \, Also note that sin 29 satisfies 


sin 2 


0-(a+-) 


cos 2 


0 -(«+-) 


sin [2 (6 — a) + 7r] 

— sin 2(6 — a), 
cos [2(6 — a) + 7r] 

— cos 2(9 — a). 


That is, the polymer orientation distribution with order parameter (cos 2(9 — a)) 
and phase angle a is identical to the polymer orientation distribution with order 
parameter — (cos 2(9 —a)) and phase angle a + n/2. So by allowing r = U (cos 2(9 — 
a)) to be both positive and negative, we can further restrict the range of a to 

r_ 7T 7T 1 

L 4 ’ 4 -I ‘ 

We first make a change of variable 9 new = 9 — a. For simplicity, we still denote 
the new variable 9 new by 9 and the corresponding pdf and potential for the new 
variable 9 new by p(9) and 1 / 1 ( 6 ), respectively. In our notation, (15) becomes 

7 * 1 ( 0 :, q, Pe) = (sin 26) = 0, (16) 

F 2 (a,q,Pe) = (cos20) = (17) 

where the new pdf p( 6 ) and the new potential ip( 6 ) are 

p( 6 ) oc exp[--0(0)] f exp [-^-9 + 1/>(9 + 6)]d9, (18) 

J 0 4 

Pe 

1 / 1 ( 6 ) = — cos 2a sin 26 — q cos 29. (19) 

For a given value of Pe, the function F\(a,q,Pe) depends on only q and a. Our 
extensive numerical calculations suggest that for any value of Pe > 0 and any value 
of a € [— f, f ], there exists one and only one value of q such that Pi (a, q , Pe) = 0. 
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A rigorous mathematical proof of the uniqueness is still open. In Appendix B, we 
prove the existence by showing the asymptotic behaviors of Pi (a, q , Pe): 


Fi(a,q,Pe) 


Pe( 1 — cos 2a)-—h • • • > 0 


Pe 

—j -I->0, 


4 q 


if cos 2a < 1, 
if cos 2a = 1, 


as q 


+00, 


Fi(a,q,Pe) = Pe(l + cos2a)-—h • • • < 0, as q 


—00. 


Since F\(a,q, Pe) is continuous and it has opposite sign at positive and negative 
infinity, there exists a point such that it is zero. In other words, equation (16) 
can always be satisfied and for a fixed value of Pe it defines uniquely a function 
q(a) such that F\{a,q{a),Pe) = 0. Notice that in (16)—(19), a appears only in the 
potential ij}( 8 ) in the form of cos 2a. Since cos 2a is an even function with respect 
to a, the whole problem (16)-(19) is even with respect to a. Consequently, q(a) is 
an even function of a. 

Now, in potential given in (19), we replace q by function q(a). It follows 
that for a fixed value of Pe, the pdf p(0) given in (18) is now a function of a only 
for a e [-f, f]. 

For mathematical convenience, let us introduce several functions as follows: 


Pe 

r(a) = q(a) + sin 2a, 

s(a) = P 2 (a,g(a),Pe), 

, x _ q(a) 
w 0 {a) = 

s{a) 

, r(a) , Pe 
w{a) = —— = Wo (a) + — 
s(a) 4 


sin 2 a 
s(a) ' 


As we discussed earlier, the function s(a) is also an even function of a. As we will 
see later, the nonlinear system (16) and (17) governing the steady state solutions 
can be conveniently described by these introduced functions. Furthermore, the non¬ 
linear system (16) and (17) in terms of these introduced functions can be simplified 
significantly. 

We plot functions q(a), s(a), wo(a) and w(a) for several values of Pe, respec¬ 
tively in Figures 1-2. 

It is clear from Figures 1-2 that function q(a), together with functions s(a) 
and wo (a) are not sensitive to the value of Pe. In contrast, the function w(a ) 
is highly dependent on Pe since the second term in the definition of w(a ) has 
a coefficient Figure 2 shows that for small values of Pe, the function w(a) 
initially increases to a local maximum; after the local maximum w(a ) decreases to 
a local minimum; after the local mininum w(a ) starts to increase again. As the 
value of Pe increases, the locations of the maximum and minimum get closer. In 
other words, as the value of Pe increases, the zigzag part of w(a ) shrinks. At a 
critical value of Pe c = 0.746027, both the value and the location of the maximum 
coincide with those of the minimum, and the zigzag part of w(a ) disappears. For 
Pe > Pe c , function w(a) is monotonically increasing for a £ [ — f • f ]• 

Note that if we replace q in equation (17) by q(a) (and replace r by r(a)), then 
equation (17) is an equation of a only where Pe and U are treated as parameters. 
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- 1/4 0 1/4 

a/n 


Figure 1. Functions q(a), s(a) and Wo(a) for various values of Pe. 




Figure 2. Function w(a) for various values of Pe. The right panel 
is a detailed view of the left panel. 

In terms of function w(a), equation (17) becomes 

w(a, Pe) = U. (20) 

Here we write it as w(a, Pe) to show explicitly its dependence on Pe. 

The phase diagram (steady states as functions of U or Pe or both) can be ob¬ 
tained by solving equation (20). As we discussed above, for Pe < Pe c , equation 
(20) has one solution for small values of U; it has three solutions for U in an in¬ 
termediate range; and it has one solution for large U. For Pe > Pe c , the zigzag 
part of function w(a) disappears. As a result, for Pe > Pe c , equation (20) has one 
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solution for any value of U >0. Figure 3 shows the phase diagram (steady state as 
a function of U) for Pe = 0.01, corresponding to the case of weak shear. As shown 
in the figure, a cusp (actually it is smooth, see the insert in Figure 3) separates the 
diagram into two branches. As we argued in [31], the branch shown in solid line, 
which corresponds to the stable branch in the case of no shear (Pe = 0), is stable. 
The branch shown in dashed line is unstable [31]. The lower part of the dashed 
line branch corresponds to the unstable isotropic branch for U > 2 in the case of 
no shear (Pe = 0). The stability of the upper part of the dashed line branch was 
analyzed in our multi-scale asymptotic study [31]. In the case of weak shear (i.e. 
Pe << 1), a nematic equilibrium will keep its shape. But under the influence of 
weak shear (Pe << 1), its phase angle a will change with respect to the slow time 
scale Ti = Pe • t. The evolution of the phase angle a is governed by the differential 
equation [31] 

da , , o , . . 

— =d(sm 2 a + c 2 ), (21) 

where C\ > 0. For — 1 < C 2 < 0, the phase angle will converge to a steady state. 
The stability of a steady state phase angle ao is determined by the derivative of 
sin 2 a + C 2 at ao- 

— (sin 2 a + c 2 ) = 2 sin a cos a = sin 2a. 
da 



Figure 3. Phase diagram for Pe = 0.01. The insert is a magnified 
view of the cusp, which separates the phase diagram into two 
branches. 

If sin2ao > 0, then ao is unstable. If sin2ao < 0, then ao is stable. As shown 
in Figure 4, for any — 1 < C 2 < 0, there are two steady state phase angles in (0,7r). 
The one in (n/2, n) is stable while the one in (0,7r/2) is unstable. The upper part of 
the dashed line branch in Figure 3 corresponds to the unstable steady state phase 
angle, and thus is unstable. In the asymptotic study of weak shear (Pe << 1) 
both the stable branch and the upper part of the unstable branch cease to exist for 
U > U c = 2.41144646. For U > U C: there is no steady state solution. Instead, there 
is a tumbling solution. In the analytical/numerical study here, the stable branch 
and the upper part of the unstable branch are connected to each other at the cusp 
(see the insert in Figure 3). 
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Figure 4. Stable and unstable steady state phase angles. 

In Figure 5 we plot the phase diagram (steady state as a function of U) for several 
values of Pe. 



Figure 5. Left panel: Phase diagram for small-to-intermediate 
values of Pe. Right panel: Phase diagram for large values of Pe. 


As shown in Figure 5, the cusp, which separates the phase diagram into two 
branches, is less and less well defined as the value of Pe increases. For Pe > Pe c , 
the cusp disappears completely. Here it is important to point out that the system 
is not driven by a potential field because the effect caused by the shear flow is not a 
potential field. As a consequnce, a steady state branch may change its stability at 
any point even when it remains isolated. In contract, when a system is driven by 
a potential field, a steady state branch may change its stability only if some other 
branch bifurcates from it (or a fold occurs) [8]. 


4. Stability study. For a given value of Pe, a steady state solution is com¬ 
pletely specified by the value of a. Once the value of a is known, U and the order 
parameter s are given by 

U = w(a), s = s(a). 

The probability density function (pdf) is calculated as follows: 


p{6) oc exp[— ip(6)\ / exp 

Jo 


-—e + ^e + e) 


de, 
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where the potential i/j(9) is 


Pe 

ip(9) = — cos 2a sin 29 — q{a) cos 29. 

Since the pdf satisfies /J r p(9)d9 = 1 and is periodic with period 7r, it has the Fourier 
series expansion of the form 


P(0) = 


M 


+ cos 2 k9 + sin2fc0) 


k =1 


where b ^ = 0. We consider perturbations of the form 


M 


Ap(9, t) = — [a*; (t) cos 2 k9 + bk (t) sin 2k9]. 

TT ‘ ^ 


( 22 ) 


(23) 


fe=l 


The Maier-Saupe interaction potential induced by the perturbed pdf p(9 , t) = p(9) + 
eAp{6,t) is 

Vms(0) = — U [(cos 29) cos 2 9 + (sin 29) sin 29 ], 

where 


cos 29) = 

f cos 29 \p{9) + eAp(9, t)\ d9 
Jo 

= 

+ eai(t), 

(sin 29) = 

f sin 29 [p{9) + eAp(9, f)] d9 
Jo 


= ebi(t). 

Substituting into the Maier-Saupe potential yields 

Vms{ 9) = —Ua < f' > cos2 9 — eU[ai(t) cos 29 + b\{t) sin 29] 

= Vms(0) + eA Vms(9), 

where Vms(9) = — Ua^ cos 29 is the unperturbed Maier-Saupe interaction po¬ 
tential and AVms[ 9) = — U[a\(t) cos20 + &i(t)sin20] is the perturbation to the 
Maier-Saupe potential. The perturbed pdf satisfies 


8 {p + eA p) 8 

dt = 89 


Y + h(0) + Vms(@) 

8 


+ £ ^Vms(®) ) ( P + £ ^p) + ■qq(p + £ Ap) 


Thus, A p(9 1 1 ) satisfies the linearized equation 


dAp 8 
= 89 


dt 

where 


TT 4 + V MS (0) ) ^P + AV MS (6)p(0) + 


dAp 


89 


, (24) 


——f — V sh (9) + V ms (6 ) — —-—I—— cos 2a cos 29 + 2q(a) sin 20,(25) 


A V M s{ 9) = 2U[a\{t) sin 20 — b\(t) cos 20]. 


(26) 
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To further evaluate the terms in (24), we appeal to (23) and arrive at 


M 


cos20-Ap = cos 20 • — [a k (t) cos2k9 + b k (t) sin2£;0] 

7r z ' 

k= 1 

2 / 1 

= — ( a/c(t)-[cos2(fc + 1)0 + cos2(fc — 1)0] 

7T z —' V 2 

k— 1 v 


+ ^fc(t) - [sin2(fc + 1)0 + sin2(fc — 1)0] 


= — [ai (t) + 02 (t) cos 20 + 6 2 (t) sin 20] 


M 


+ -V^ ([afe_i(t) + a k +i(t)} cos2fc0 + [6fe_i(t) + 6 fc+ i(t)] sin2A:0). 

7T z ' 


k—2 


Similarly, we have 


sin20-Ap = — [bi(t) + 6 2 (t) cos 20 — 02 (t) sin 20] 


M 


H— ([0fc+i(t) - 0fc-i(t)] cos2 k6 + [a k -i(t) - a k +i(t)] sin2fc0). 

7T ^^ 


k—2 


Using the expansion (22) and noting that b^ r> = 0, we find after some straightfor¬ 
ward manipulations that 


sin 20 ■ p = sin 20 • 


M 


+ y^(t4° cos 2 kd + sin2fc0) 


fc =1 


1 r 


cos 20 + (1 — a'' 2 > ) sin 20 


(Oh 


1 M 

~ [( 6 fc+i “ b k-i) cos 2fc0 + (4°^ - 4°i) sin2fc0 


k=2 


and 


1 r 


cos 20 -p = — 4°^ + (1 + cos 20 + b^ 1 sin 20 


(0)x 


1.(0), 


1 M 

- 5Z [( a i-i+4+i) cos 2ke +(4°-i+4+i) sin 2ke 


k—2 
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From (22), 



I = 


and 

AVms(0)p(0) 


, (23), , (25) 

Vsh(0) + Vms(&) 


and 
A p as 


(26), 


we express 


Pe 

T 


Pe 


cos 2 a cos 20 + 2 q(a) sin 26 


A p 


1 

7T 


Pe 

T 


Pe 


cos 2a ■ a i (i) + 2 g(a) • £>i (t) 


1 



Pe 

-Pe ■ a± (t) + — cos 2a ■ a 2 (i) + 2q(a) • b 2 (i) 


cos 29 


1 



Pe 

-Pe-bi (t) + — cos 2a ■ b 2 (t) — 2q(a) ■ a 2 (t) 


sin 29 


1 M / p e 

H— V -Pe • a fe (t) + — cos 2a ■ [a k -i(t) + a k+1 (t )] 

7T z ' \ Z 

k=2 v 


+ 2 <?(a) • [ 6 fc+ i(t) - 6 fc_i(t)]) cos 2 fc 0 

1 M / p ^ 

H— V ( -Pe • bk(t) + —^ cos 2 a • [b k -i{t) + b k+1 (t)\ 

7T z ' \ Z 

/c—2 v 

+ 2 g(a) • [a k -i(t) - a fe+ i(t)]) sin 2 fc 0 , 


= 2U[a\(t) sin20 — b\(t) cos20]p(0) 


2U_ 

7T L 

2 U 

7T 

2U 


M 


+ 


+ - 


7 r 


-fei(i) • ai 
ai(t) • - 61 (i) • (1 + a 2 UJ ) 

ai(t) ■ (1 - a' 0) ) - b^t) ■ & 2 ° } 


cos 20 

sin 20 


M 


+ V 51 [ ai (*) • “ b k- 1 ) - ’ ( a i-i + a i+i) 

k—2 
M 


+v £ M*) ■ ( a i-i~ a kh) - • ( & i°-i+ 6 L+i) 


k—2 


Pe 

qi = — cos 20 , g 2 = 2 q(a), g 3 = 2C7. 


quantity 


cos 2 fc 0 


sin 2 fc 0 . 


Let 
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Putting all these in the right hand side of (24) and simplifying, we get 
2 r 

RHS = - Pe ■ bi(t) + qib 2 {t) - q 2 a 2 {t) - q 3 [a { 2 ^ - 4, 0) ]ai(t) 


-'73[&o° ) + b 2° ) ] i »l( t ) - 4oi(t) 


cos 20 


+ - 


2 r 


Pe ■ ai(t ) - qia 2 (t) - q 2 b 2 {t) - q 3 [b ( 2 0) - 4°Vi( i ) 


+93 [4° } + 4 0) ] & 1 W - 4& i(i) 


x 0 

M 


sin 2 6 


H— V' fc [—-Pe • 6 fe (t) + qi[b k -i(t) + b k+1 (t)} - q 2 [a k+1 (t) - a k -i{t)\ 

7r ^^ 


k=2 

,(°) „(0) 


-93[afc+! - 4-i] ’ «i(*) “ 93 [4- 1 + 4+i] ' & iW - 4fca fcW 


cos 2/c# 


M 


H— V' k [Pe ■ a k (t) - qi[a k -i{t) + a k+1 (t)[ - q 2 [b k+1 (t ) - & fe _i(t)] 

7T z ^ 


k—2 


Hence, from 
da\{t) 
dt 


dbi(t) 

dt 


da k (t) 

dt 


db k (t ) 
dt 


-93[4+i - 4-J • ai(i) + 93[a£-i + 4+iI ’ 

—4/cfofc(t)] sin 2fc(9. (27) 

(24), it follows that the Fourier coefficients of the perturbation satisfy 

= -Pe • foi(i) + q\b 2 {t) - q 2 a 2 {t) - q^a^ - a^ViW 
-93 [&o 0) + £ > 2° ) ] £> i (t) - 4ai (t), 

= Pe • ai(i) - gia 2 (t) - q 2 b 2 (t) - q^b^ - 4 0) ]ai(i) 

+93[ao 0) + 4 0) ] & i(0 - 4*i(*)> 

= k [-Pe ■ b k (t ) + qi[b k _i(t) + b k+1 (t)] - q 2 [a k+1 (t) - a k -i(t )] 

—73[a£+i - 4-iM*) - 93 [4-i + 4+ilM*) - 4ka k {t) , 

= k [Pe ■ a k {t) - 9 i[a fc _i(t) + a k+ i(t)] - q 2 [b k+1 (t) - & fc _i(t)] 

-93[4+1 - 4-iM*) + 93(4-! + a kh]bi(t) - 4 kb k (t)] , 


where we have set 


4 0) = i, 4 0) = o. 


If we collect all of the Fourier coefficients of A p(6, t) into a vector 
u (t) = [a 1 (t),b 1 {t),a 2 {t),b 2 (t),-■ • , a M {t), b M {t)] T , 
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then the above equations can be cast in a compact form 

(28) 

where we have used the notation A[p{6)\ to show explicitly that the matrix A 
depends on the unperturbed steady state p(9). If all of the eigenvalues of A have 
negative real parts, then the perturbation decays to zero as time goes to infinity 
and the steady state p{6) is linearly stable; if at least one of the eigenvalues of A 
has a positive real part, then the steady state p(6) is unstable. We use the Matlab 
package to solve this eigenvalue problem to determine the stability. 

5. Numerical results. Referring to Figure 3, we showed the phase diagram 
with stability/instability indicated for the case of Pe = 0.01. In the figure, the solid 
line represents the stable branch and the dashed line represents the unstable branch. 
The numerical stability result is consistent with that of the asymptotic analysis. 
Numerical stability analysis shows that the stable and the unstable branches are 
connected at the folding point of the phase curve where the phase curve is vertical 
(i.e., the point where ds/dU = oo). The transition point from the stable branch to 
the unstable branch can only be calculated from the numerical stability analysis. 
The asymptotic analysis is not capable of resolving this detail. Figure 6 shows the 
phase diagram with stability/instability indicated for the case of Pe = 0.1. The 
solid line represents the stable branch and the dashed line represents the unstable 
branch. Numerical stability analysis shows that the transition point from the stable 
branch to the unstable branch is the folding point of the phase curve where the phase 
curve is vertical (i.e. the point where ds/dU = oo). 



Figure 6. Phase diagram for Pe = 0.1. 

In Figure 7 we show the phase diagrams with stability/instability indicated for 
several intermediate values of Pe. The solid line represents the stable branch and 
the dashed line represents the unstable branch. Numerical stability analysis shows 
that the transition from the stable branch to the unstable branch is the folding 
point of the phase curve if such a folding point exists. Note that this qualitative 
behavior based on the closure models has been observed in [21] even though they 
differ qualitatively. 
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Figure 7. Phase diagram for intermediate values of Pe. 

For Pe < 0.746, such a folding point exists. For Pe > 0.746, there is no such 
folding point. For Pe > 0.746, the transition from the stable branch to the unstable 
branch occurs at a point beyond the maximum. 



Figure 8. Phase diagram for large values of Pe. 


Figure 8 depicts the phase diagram with stability/instability indicated for sev¬ 
eral large values of Pe. The solid line corresponds to the stable branch whereas 
the dashed line corresponds to the unstable branch. Numerical stability analysis 
indicates that the transition from the stable branch to the unstable branch occurs 
at a point beyond the maximum. Furthermore, as the value of Pe increases, the 
value (U) of the transition point increases. 

Figure 9 describes the region of stable steady state (shaded) and the region of 
tumbling solution in the (U, Pe)-plane. This figure is also inherent in [21], where 
the locus of the Hopf bifurcation was shown. Here, it is captured from the kinetic 
theory. 
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Figure 9. Region of stable steady state solution (shaded) and the 
region of tumbling solution. Left panel: overall view; Right panel: 
detailed view for small Pe. 


6. Conclusions. We have semi-analytically investigated all steady state solu¬ 
tions and their stability for planar 2D nematic polymers under an imposed shear 
flow with arbitrary strength. We have confirmed the mesoscopic closure models of 
Lee et al. [21], extending their analysis and numerical studies to the Smoluchowski 
equation for the orientational probability distribution. A detailed phase diagram is 
given, consisting of all stable steady states versus the normalized concentration U 
and normalized shear rate or Peclet number Pe. Special attention has been given to 
the locus of bifurcations in the phase diagram where stable steady states no longer 
exist. These bifurcations are shown to consist of a fold among stable and unstable 
branches up to some critical Pe, after which the fold disappears and is replaced by 
Hopf bifurcation. Extension of the current work to 3-D model will be much more 
challenging mathematically. 

Appendix A An efficient method for calculating f n exp [— fd + ijj(8 + (9)] dd 
Using the Fourier series expansion of exp(^(0)) = YJJk=-oo c fc exp(ifc20), we ob¬ 
tain 


p7T °° /*7T 

/ exp [— fd + ip(d + 9)] dd = V Cfc / exp [— fd + ik2(d + 9)] dO 

Jo „ Jo 


where 


Uo 


exp [(—/ + ik2)0\ dd 


k= — oo 
oo 

Ck 

k=— oo 




exp(ik29) 


k— — oo 


/ — ik2 


dk — Cfc 


dk exp(*A:20), 

k— — oo 

1 - exp(—/ t) 


(29) 


/ — ik2 

The numerical procedure for calculating exp [—fd + tjj(8 + 0)] dd is as follows. 
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Apply FFT on ip(9) to calculate c k . 

Calculate <4 = c k -—' 

Apply IFFT on d k to compute exp \—f0 + 1 (>(9 + 0)] cL9. 


Appendix B Existence of solution of F\(a,q, Pe) = (sin2 9) = 0 
We consider several cases. 

Case 1: q —> +oo 

Since the pdf p(9) given in (18) is periodic with period ir , we can consider any 
interval of length n. Without loss of generality, let us select the interval [—§,§]■ 
We have 


Fi(a,q,Pe) = (sin 29) 


where 


f-ir 


Z , 

1 


sin 2 9 exp 




--9 + ^(9 + 9)-m 


d9d9 


f ■Vr i-TT 


sin 2 9 exp 


- 5 - 


Pc — Pc — 

——9 + — cos2a[sin2(0 + 9) — sin 29] 


— g[cos 2(9 + 9) — cos 29]\ d9d9, 


(30) 


Z 



Pc — Pc — 

——6 + — cos 2a[sin 2(9 + 9) 


— gfcos 2(9 + 9) — cos 20]] d9d9. 


sin 20] 


For positive and large q , the dominant contribution in the above two integrals 
comes from a small region near 0 = 0 and 0 = ^. Away from this small region, the 
contribution is exponentially small. 

In the inner integral with respect to 0, we make a change of variables 9 new = 0+0. 
For simplicity, we still use 0 to denote the new variable 9 new . Since the dominant 
contribution comes from a small region near 0 = 0 and 0 = ((, we can extend the 
integration limits away from the region of the dominant contribution. Specifically, 
we rewrite (30) as 
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Fi(a,q,Pe) 

1 Pe 

— / sin20exp[—0 


Pe 

cos 2a sin 20 + g cos 20] 


x 



— ^(0 + 9) + ^ cos2asin2(0 + 9) 


geos 2(0 + 9) 


d9d9 


1 

~Z 


Pc Pc 

sin 29 exp[-^-0-— cos 2a sin 2 9 + q cos 20] 


x 



Pe - Pe - - 

/ exp 

Je 

——9 4—— cos 2a sin 29 — q cos 20 
2 4 


d9d9 + T.S.T. 


1 f* Pe 

— / sin20exp[—0 


Pe 

cos 2a sin 29 + q cos 20] 


x / exp 

J o 


Pc - Pc - - 

——9 -\—— cos 2a sin 20 — q cos 20 
2 4 


d9d9 


+T.S.T. 

Here T.S.T. stands for transcendentally small term. Let 
Ci = 


We have 


/*71 

Jo 


Ci 


exp 


Pc - Pc - - 

-0 H-cos 2a sin 29 — q cos 20 

2 4 


d6. 


(31) 

(32) 


_ Pc Pc 

Fi(a,q, Pe) = ^- / sin20exp [— 0 — cos 2a sin 20 + geos 20] + T.S.T. 


Z _ 

2 Ci 

~Z~ 
+T.S.T. 


' Pe Pe 

sin 20 sinh ( —0-cos 2a sin 20 1 exp(g cos 20) 


(33) 

We apply the Laplace method to the integral in (33). The dominant contribution 
comes from a small region near 0 = 0. To proceed, we need to find the leading term 
expansion of the part not containing q around 0 = 0. We discuss two cases. 

Case 1A : cos 2a < 1 
In this case, we have 

f Pe Pe \ 

sin 20 sinh ( -^-9 -— cos 2a sin 20 J = Pe(l — cos 2a)0 2 + • • • , 


and 


Fi(a,q,Pe) 



Pe( 1 — cos 2a)0 2 exp 



MV 

2 ' 


d9 + • • • 


= exp(g)Pe(l 


cos 2a) 



0 2 exp[— 2q9 2 ]d9 -\ - 
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Introducing a change of variables 

, 1 1 
s = 6 , dO = — ds = —=ds, 

29 2y/s 

we find 

C f 00 

Fi(a,g, Pe) = —- exp(g)Pe(l — cos 2a) / v^ ex P[ — H - 

% J o 

= Tjr exp(g)Pe(l - cos 2a) + • • • ■ (34) 

Similarly, we get 

Z = 2Ci f cosh(^-0 — —p cos2asin20) exp(gcos20)d0 + T.S.T. 

Jo 2 4 

/»00 

= 2Ci exp(q) / exp(— 2q9 2 )d8 + • • • 

Jo 


= Ciexp(q)-^L H-. (35) 

Substituting (35) into (34) yields 

Fi(a, g, Pe) = Pe(l — cos 2a)— ->0, as q —> +oo. 

4g 

Case IB : cos 2a = 1 

When cos 2a = 1, we have 

P & Pe ^Pe 

sin2#sinh(-^-$-— cos 2 a sin 29) = —0 4 + • • • 

and it follows that 

Fi(a,q,Pe) = 


Combining (35) and (36), we obtain 

Pe 

Fi(a,q, Pe) = —- H-> 0, as q —> +oo. 

X/7 Z 


9 Pp r°° 

^T-exp(g)— J 9 4 exp(—2q9 2 )dd + ■ 


Ci . , 2Pe f°° . n , , 

■y e xp(g)-g- J svsexp(—2gs)as 

Cr , ,2Pe |V5F 
Z 6XP(9) 3 (2q) 2 ^/2q " 


(36) 


Case 2 : q —> —oo 

This case is similar to Case 1 but simpler than it. The derivation is skipped. For 
negatively large q, the asymptotic result is 


Pi (a, q, Pe) = Pe(l + cos 2a) ——h ■ 


<0, as q —> —oo, 


where we use the fact that the coefficient (1 + cos 2a) is always positive for a £ 

7T 7r 

4 ’ 4 
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